#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;
static doublereal c_b15 = 10.;
static doublereal c_b88 = 0.;

/*   VERSION 2    DOES NOT USE EISPACK */

/* ------------------------------------------------------------------ */

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::dnlaso_(U_fp
   op, U_fp iovect, integer *n, integer *nval, integer *nfig, integer *nperm,
   integer *nmval, doublereal *val, integer *nmvec, doublereal *vec, integer
   *nblock, integer *maxop, integer *maxj, doublereal *work, integer *ind,
   integer *ierr)
{
  /* System generated locals */
  integer vec_dim1, vec_offset, val_dim1, val_offset, i__1, i__2;
  doublereal d__1;

  /* Local variables */
  static doublereal temp, tarr[1];
  //extern doublereal dnrm2_(integer *, doublereal *, integer *);
  static integer i__, m, nband;
  static doublereal delta;
  //extern /* Subroutine */ int dnwla_(U_fp, U_fp, integer *, integer *, 
  // integer *, integer *, integer *, doublereal *, integer *, 
  // doublereal *, integer *, integer *, integer *, integer *, 
  // doublereal *, doublereal *, doublereal *, doublereal *, 
  // doublereal *, doublereal *, doublereal *, doublereal *, 
  // doublereal *, doublereal *, doublereal *, doublereal *, 
  // doublereal *, doublereal *, integer *, logical *, logical *, 
  // doublereal *, doublereal *, integer *);
  static logical small;
  //extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
  // doublereal *, integer *);
  static integer i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, nv;
  //extern /* Subroutine */ int dnppla_(U_fp, U_fp, integer *, integer *, 
  // integer *, integer *, doublereal *, integer *, doublereal *, 
  // integer *, doublereal *, doublereal *, doublereal *, doublereal *,
  //  doublereal *, doublereal *, doublereal *, logical *, logical *, 
  // doublereal *);
  static logical raritz;
  //extern /* Subroutine */ int dortqr_(integer *, integer *, integer *, 
  // doublereal *, doublereal *), dvsort_(integer *, doublereal *, 
  // doublereal *, integer *, doublereal *, integer *, integer *, 
  // doublereal *);
  static doublereal eps;
  static integer nop;

  /* AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT */
  /* COMPUTER SCIENCES DEPARTMENT */
  /* UNIVERSITY OF TEXAS AT AUSTIN */
  /* AUSTIN, TX 78712 */
  /* VERSION 2 ORIGINATED APRIL 1982 */
  /* CURRENT VERSION  JUNE 1983 */
  /* DNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF */
  /* THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX.  THE SUBROUTINE */
  /* DNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE DNWLA WHICH IMPLEMENTS */
  /* THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND */
  /* SUBROUTINE DNPPLA WHICH POST PROCESSES THE OUTPUT OF DNWLA. */
  /* HOWEVER DNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING */
  /* PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS. */
  /* DNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES.  IF */
  /* THE LARGEST EIGENVALUES ARE DESIRED DNLASO IMPLICITLY USES THE */
  /* NEGATIVE OF THE MATRIX. */

  /* ON INPUT */


  /*   OP   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE */
  /*     OP(N,M,P,Q).  P AND Q ARE N X M MATRICES AND Q IS */
  /*     RETURNED AS THE MATRIX TIMES P. */

  /*   IOVECT   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE */
  /*     IOVECT(N,M,Q,J,K).  Q IS AN N X M MATRIX.  IF K = 0 */
  /*     THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH */
  /*     THE JTH LANCZOS VECTORS.  IF K = 1 THEN Q IS RETURNED */
  /*     AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS.  SEE */
  /*     DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES. */

  /*   N   THE ORDER OF THE MATRIX. */

  /*   NVAL   NVAL SPECIFIES THE EIGENVALUES TO BE FOUND. */
  /*     DABS(NVAL)  IS THE NUMBER OF EIGENVALUES DESIRED. */
  /*     IF NVAL < 0 THE ALGEBRAICALLY SMALLEST (LEFTMOST) */
  /*     EIGENVALUES ARE FOUND.  IF NVAL > 0 THE ALGEBRAICALLY */
  /*     LARGEST (RIGHTMOST) EIGENVALUES ARE FOUND.  NVAL MUST NOT */
  /*     BE ZERO.  DABS(NVAL) MUST BE LESS THAN  MAXJ/2. */

  /*   NFIG   THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE */
  /*     EIGENVALUES.  NFIG MUST BE GREATER THAN OR EQUAL TO 1. */

  /*   NPERM   AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER */
  /*     SUPPLIED EIGENPAIRS.  IN MOST CASES NPERM WILL BE ZERO.  SEE */
  /*     DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER */
  /*     THAN ZERO.  NPERM MUST NOT BE LESS THAN ZERO. */

  /*   NMVAL   THE ROW DIMENSION OF THE ARRAY VAL.  NMVAL MUST BE GREATER */
  /*     THAN OR EQUAL TO DABS(NVAL). */

  /*   VAL   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW */
  /*     DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4.  IF NPERM */
  /*     IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED */
  /*     IN VAL.  SEE DOCUMENTATION FOR DETAILS. */

  /*   NMVEC   THE ROW DIMENSION OF THE ARRAY VEC.  NMVEC MUST BE GREATER */
  /*     THAN OR EQUAL TO N. */

  /*   VEC   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW */
  /*     DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST DABS(NVAL).  IF */
  /*     NPERM > 0 THEN THE FIRST NPERM COLUMNS OF VEC MUST */
  /*     CONTAIN THE USER SUPPLIED EIGENVECTORS. */

  /*   NBLOCK   THE BLOCK SIZE.  SEE DOCUMENTATION FOR CHOOSING */
  /*     AN APPROPRIATE VALUE FOR NBLOCK.  NBLOCK MUST BE GREATER */
  /*     THAN ZERO AND LESS THAN  MAXJ/6. */

  /*   MAXOP   AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE */
  /*     OP.  DNLASO TERMINATES WHEN MAXOP IS EXCEEDED.  SEE */
  /*     DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP. */

  /*   MAXJ   AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND */
  /*     DOCUMENTATION ON IOVECT).  FOR THE FASTEST CONVERGENCE MAXJ */
  /*     SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE */
  /*     MAXJ LARGER THAN MAXOP*NBLOCK. */

  /*   WORK   A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS */
  /*     LARGE AS */

  /*         2*N*NBLOCK + MAXJ*(NBLOCK+NV+2) + 2*NBLOCK*NBLOCK + 3*NV */

  /*            + THE MAXIMUM OF */
  /*                 N*NBLOCK */
  /*                   AND */
  /*         MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1) */

  /*     WHERE NV = DABS(NVAL) */

  /*     THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED */
  /*     STARTING VECTORS.  SEE DOCUMENTATION FOR GUIDELINES IN */
  /*     CHOOSING STARTING VECTORS. */

  /*   IND   AN INTEGER ARRAY OF DIMENSION AT LEAST DABS(NVAL). */

  /*   IERR   AN INTEGER VARIABLE. */


  /* ON OUTPUT */


  /*   NPERM   THE NUMBER OF EIGENPAIRS NOW KNOWN. */

  /*   VEC   THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS. */

  /*   VAL   THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING */
  /*     EIGENVALUES.  THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF */
  /*     THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN- */
  /*     VALUES.  THE THIRD COLUMN CONTAINS MORE DOUBLE PRECISIONISTIC ESTIMATES */
  /*     OF THE ACCURACY OF THE EIGENVALUES.  THE FOURTH COLUMN CONTAINS */
  /*     ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS.  SEE */
  /*     DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES. */

  /*   WORK   IF WORK IS TERMINATED BEFORE COMPLETION (IERR = -2) */
  /*     THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS */
  /*     FOR RESTARTING THE ALGORITHM AND DNLASO CAN BE IMMEDIATELY */
  /*     RECALLED TO CONTINUE WORKING ON THE PROBLEM. */

  /*   IND   IND(1)  CONTAINS THE ACTUAL NUMBER OF CALLS TO OP.  ON SOME */
  /*     OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER */
  /*     THAN MAXOP. */

  /*   IERR   AN ERROR COMPLETION CODE.  THE NORMAL COMPLETION CODE IS */
  /*     ZERO.  SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO */
  /*     COMPLETION CODES. */


  /* INTERNAL VARIABLES. */



  /* NOP   RETURNED FROM DNWLA AS THE NUMBER OF CALLS TO THE */
  /*   SUBROUTINE OP. */

  /* NV   SET EQUAL TO DABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED, */
  /*   AND PASSED TO DNWLA. */

  /* SMALL   SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED. */

  /* RARITZ   RETURNED FROM DNWLA AND PASSED TO DNPPLA.  RARITZ IS .TRUE. */
  /*   IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED. */

  /* DELTA   RETURNED FROM DNWLA AS THE EIGENVALUE OF THE MATRIX */
  /*   WHICH IS CLOSEST TO THE DESIRED EIGENVALUES. */

  /* DNPPLA   A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED */
  /*   BY DNWLA. */

  /* DNWLA   A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM */
  /*   WITH SELECTIVE ORTHOGONALIZATION. */

  /* DMVPC   A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND */
  /*   ORTHOGONALITY COEFFICIENT OF GIVEN RITZ VECTORS. */

  /* DORTQR   A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS */
  /*   USING HOUSEHOLDER REFLECTIONS. */

  /* DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP   A SUBSET OF THE BASIC LINEAR */
  /*   ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION. */

  /* DLARAN   A SUBROUTINE TO GENERATE RANDOM VECTORS */

  /* DLAEIG, DLAGER, DLABCM, DLABFC   SUBROUTINES FOR BAND EIGENVALUE */
  /*   CALCULATIONS. */

  /* ------------------------------------------------------------------ */

  /* THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS. */

  /* Parameter adjustments */
  val_dim1 =  *nmval;
  val_offset = 1+val_dim1 * 1;
  val -= val_offset;
  vec_dim1 =  *nmvec;
  vec_offset = 1+vec_dim1 * 1;
  vec -= vec_offset;
  --work;
  --ind;

  /* Function Body */
  nv = abs(*nval);
  ind[1] = 0;
  *ierr = 0;
  if (*n <  *nblock * 6)
  {
    *ierr = 1;
  }
  if (*nfig <= 0)
  {
    *ierr += 2;
  }
  if (*nmvec <  *n)
  {
    *ierr += 4;
  }
  if (*nperm < 0)
  {
    *ierr += 8;
  }
  if (*maxj <  *nblock * 6)
  {
    *ierr += 16;
  }
  if (nv < max(1,  *nperm))
  {
    *ierr += 32;
  }
  if (nv >  *nmval)
  {
    *ierr += 64;
  }
  if (nv >  *maxop)
  {
    *ierr += 128;
  }
  if (nv >=  *maxj / 2)
  {
    *ierr += 256;
  }
  if (*nblock < 1)
  {
    *ierr += 512;
  }
  if (*ierr != 0)
  {
    return 0;
  }

  small =  *nval < 0;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS. */
  /* IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION */
  /* OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO  -1 */
  /* AND DNLASO TERMINATES. */

  if (*nperm == 0)
  {
    goto L110;
  }

  /* THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST */
  /* EIGENVALUES ARE DESIRED, SINCE DNWLA WILL IMPLICITLY USE THE */
  /* NEGATIVE OF THE MATRIX. */

  if (small)
  {
    goto L20;
  }
  i__1 =  *nperm;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    val[i__ + val_dim1] =  - val[i__ + val_dim1];
    /* L10: */
  }

  /* THIS SORTS THE USER SUPPLIED VALUES AND VECTORS. */

  L20: dvsort_(nperm, &val[val_offset], &val[(val_dim1 << 1) + 1], &c__0, tarr,
    nmvec, n, &vec[vec_offset]);

  /* THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON. */
  /* IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE. */

  i__1 =  *nperm;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    val[i__ + (val_dim1 << 1)] = (d__1 = val[i__ + (val_dim1 << 1)], abs(d__1));
    val[i__ + val_dim1 * 3] = dnrm2_(n, &vec[i__ *vec_dim1 + 1], &c__1);
    /* L60: */
  }

  /* THIS PERFORMS THE ORTHONORMALIZATION. */

  m =  *n **nblock + 1;
  dortqr_(nmvec, n, nperm, &vec[vec_offset], &work[m]);
  m =  *n **nblock -  *nperm;
  i__1 =  *nperm;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    m = m +  *nperm + 1;
    if ((d__1 = work[m], abs(d__1)) > val[i__ + val_dim1 * 3]*(float).9)
    {
      goto L70;
    }
    *ierr =  - 1;
    return 0;

    L70: ;
  }

  /* THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN */
  /* THE ARRAY WORK FOR LATER REFERENCE IN DNWLA. */

  m = (*n << 1) **nblock + 1;
  dcopy_(nperm, &val[(val_dim1 << 1) + 1], &c__1, &work[m], &c__1);

  /* THIS SETS EPS TO AN APPROXIMATION OF THE RELATIVE MACHINE */
  /* PRECISION */

  /* ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT */
  /* ***IN A PRODUCTION CODE */

  L110: eps = 1.;
  for (i__ = 1; i__ <= 1000; ++i__)
  {
    eps *= .5;
    temp = eps + 1.;
    if (temp == 1.)
    {
      goto L130;
    }
    /* L120: */
  }

  /* ------------------------------------------------------------------ */

  /* THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM */
  /* WITH SELECTIVE ORTHOGONALIZATION. */

  L130: nband =  *nblock + 1;
  i1 =  *n **nblock + 1;
  i2 = i1 +  *n **nblock;
  i3 = i2 + nv;
  i4 = i3 + nv;
  i5 = i4 + nv;
  i6 = i5 +  *maxj * nband;
  i7 = i6 +  *nblock **nblock;
  i8 = i7 +  *nblock **nblock;
  i9 = i8 +  *maxj *(nv + 1);
  i10 = i9 +  *nblock;
  i11 = i10 + (nv << 1) + 6;
  i12 = i11 +  *maxj *((*nblock << 1) + 1);
  i13 = i12 +  *maxj;
  dnwla_((U_fp)op, (U_fp)iovect, n, &nband, &nv, nfig, nperm, &val[val_offset],
    nmvec, &vec[vec_offset], nblock, maxop, maxj, &nop, &work[1], &work[i1],
    &work[i2], &work[i3], &work[i4], &work[i5], &work[i6], &work[i7], &work[i8],
    &work[i9], &work[i10], &work[i11],  &work[i12], &work[i13], &ind[1], &small,
    &raritz, &delta, &eps, ierr);

  /* ------------------------------------------------------------------ */

  /* THIS SECTION CALLS DNPPLA (THE POST PROCESSOR). */

  if (*nperm == 0)
  {
    goto L140;
  }
  i1 =  *n **nblock + 1;
  i2 = i1 +  *nperm **nperm;
  i3 = i2 +  *nperm **nperm;
  /* Computing MAX */
  i__1 =  *n **nblock, i__2 = (*nperm << 1) **nperm;
  i4 = i3 + max(i__1, i__2);
  i5 = i4 +  *n **nblock;
  i6 = i5 + (*nperm << 1) + 4;
  dnppla_((U_fp)op, (U_fp)iovect, n, nperm, &nop, nmval, &val[val_offset],
    nmvec, &vec[vec_offset], nblock, &work[i1], &work[i2], &work[i3], 
    &work[i4], &work[i5], &work[i6], &delta, &small, &raritz, &eps);

  L140: ind[1] = nop;
  return 0;
} /* dnlaso_ */
